clear 
est clear

********************************************************************************
*	Programs
******************************************************************************** 

*** program for the OLS regressions... 
capture program drop regs
program define regs, eclass

	*** Syntax 
	syntax varlist(numeric) [, treatment(varlist numeric) rounds(numlist) ctrls(varlist numeric)]
 
	*** loop over outcomes and rounds
	qui foreach y of local varlist {
		foreach r in `rounds' {
		    
			*** create the block variables, always centred for sample with nonmissing outcome 
			capture drop blocks*
			qui tab block, gen(blocks)
			qui foreach v of varlist blocks* {
				replace `v' = . 			if missing(`y') & round ==`r'
				su `v' 						if round == `r' & !missing(`y')
				replace `v' = `v' - r(mean) if round == `r' & !missing(`y')
			}
			
			*** Model, includes: 
			*	- treatment,
			*	- block FE fully interacted,
			*	- clustered SE's at club level,
			*	- tracking weights (=1 in RR and 12m fups)
			* 	- controls fully interacted (when adjusted)
			local model_name = "`y'_`r'"
			if (length("`ctrls'")>0) local model_name = "`y'_`r'_adj"
			
			eststo `model_name': reg `y' b0.`treatment'##c.(`ctrls' blocks*) [pw=w] if round==`r', vce(cluster club)
			
			*** save stats 
			* control mean and SD 
			qui su `y' [aw=w] if `treatment'==0 & round==`r' 
			estadd scalar c_mean = r(mean) 	
			estadd scalar c_sd = r(sd) 
			* test coeffs p-values, when treatments are combined
			if ("`treatment'"=="treat_ther") { 
				* p-value for IPT-G
				test 1.treat_ther == 0 
				estadd scalar pval_ther_any = r(p)
			}
			* p-value for the two treatments and comparison stats 
			else if ("`treatment'"=="treat") {
			    * p-value for IPT-G
				test 1.treat == 0 
				estadd scalar pval_ther = r(p)
				* p-value for IPT-G+
				test 2.treat == 0 
				estadd scalar pval_cash = r(p) 
				* diff between treaments 
				lincom 2.treat - 1.treat
				estadd scalar est_diff = r(estimate) 
				estadd scalar pval_diff = r(p) 
			}
			if (length("`ctrls'")==0) 		estadd local ctrls "No"
			else if (length("`ctrls'")>0)	estadd local ctrls "Yes"
		}
	}
	drop blocks* 
end


*** program to extract q values for the outcomes and treatment specified, within each round
capture program drop get_qvals
program define get_qvals, eclass

	*** Syntax 
	syntax varlist(numeric) [, rounds(numlist) tname_in(string) tsuff_out(string) method(string) adj(str)]
	
	*** Create vector of p-values
	* matrix to store them
	local varcount : word count `varlist' 
	matrix pvals_tab = J(`varcount',1,.)  
	matrix colnames pvals_tab = pvals
	
	* loop over outcomes and rounds
	qui foreach r in `rounds' {
	    local i = 0
		foreach y of local varlist {
		    local ++i
			* restore estimates
			est restore `y'_`r'`adj'
			est replay
			* collect pvals in matrix 
			if ("`tname_in'" == "est_diff") {
				lincom 2.treat - 1.treat 
				matrix pvals_tab[`i',1] = `r(p)'
				matrix list pvals_tab				
			}
			else {
				matrix list r(table)
				matrix pvals_tab[`i',1] = r(table)["pvalue","`tname_in'"]
				matrix list pvals_tab			    
			}
		}
		* save in a dataset 
		capture drop pvals`r'
		svmat pvals_tab, names(col)
		rename pvals pvals`r'
		
		*** get q-values  
		capture drop qvals`r'
		qqvalue pvals`r', method(`method') qvalue(qvals`r') 
	}
	
	*** add to estimates 
	* loop over outcomes and rounds
	qui foreach r in `rounds' {
	    local i = 0
		foreach y of local varlist {
		    local ++i
			* restore estimates
			est restore `y'_`r'`adj'
			* add as a scalar
			local qv = qvals`r' in `i'
			di `qv'
			estadd scalar qval_`tsuff_out' = `qv'
		}
	} 
end




********************************************************************************
*	I - Unadjusted results in paper
********************************************************************************

use "${data}/Analysis_wide.dta", clear
distinct block club

*** prep the data 
* get it in long format
reshape long phq8_score w, i(block cr_id club treat ) j(round)
* keep only what we need 
keep block cr_id club treat treat_ther round w ///
	 age_b everpregnant_b evermarried_b PPI_score_b phq8_score_b ///
	 phq8_score 
* drop missings, use the phq score that was always asked of everyone
drop if missing(phq8_score)

*** create alternative indicators 
gen phq8_mdd = phq8_score >= 10 if !missing(phq8_score)
gen phq8_severe = phq8_score >= 15 if !missing(phq8_score)


*** Regressions
* RR
regs phq8_mdd phq8_severe, ///
	treatment(treat_ther) rounds(1) 
regs phq8_mdd phq8_severe, ///
	treatment(treat_ther) rounds(1) ctrls(phq8_score_b age_b everpregnant_b evermarried_b PPI_score_b)
	
	
* 12 and 24 months
regs phq8_mdd phq8_severe, ///
	treatment(treat) rounds(2 3) 
regs phq8_mdd phq8_severe, /// 
	treatment(treat) rounds(2 3) ctrls(phq8_score_b age_b everpregnant_b evermarried_b PPI_score_b)


*** MHT adjustments  
*unadjusted 
get_qvals phq8_mdd phq8_severe, rounds(1) tname_in("1.treat_ther") tsuff_out("ther_any") method("simes")
get_qvals phq8_mdd phq8_severe, rounds(2 3) tname_in("1.treat") tsuff_out("ther") method("simes")
get_qvals phq8_mdd phq8_severe, rounds(2 3) tname_in("2.treat") tsuff_out("cash") method("simes")
get_qvals phq8_mdd phq8_severe, rounds(2 3) tname_in("est_diff") tsuff_out("diff") method("simes")
*adjusted 
get_qvals phq8_mdd phq8_severe, rounds(1) tname_in("1.treat_ther") tsuff_out("ther_any") method("simes") adj("_adj")
get_qvals phq8_mdd phq8_severe, rounds(2 3) tname_in("1.treat") tsuff_out("ther") method("simes") adj("_adj")
get_qvals phq8_mdd phq8_severe, rounds(2 3) tname_in("2.treat") tsuff_out("cash") method("simes") adj("_adj")
get_qvals phq8_mdd phq8_severe, rounds(2 3) tname_in("est_diff") tsuff_out("diff") method("simes") adj("_adj")


	
*** Table with adjusted and unadjusted
* RR:
#delimit ;
esttab 	phq8_mdd_1 phq8_mdd_1_adj phq8_severe_1 phq8_severe_1_adj
	using "${tables}/results-mh-mdd-revised-1.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat_ther) 
	stats(c_mean c_sd N ctrls pval_ther_any qval_ther_any , fmt(3 3 0 0 3)  
	label("Control mean" "Control SD" "Observations" "Controls" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* 12m:
#delimit ;
esttab 	phq8_mdd_2 phq8_mdd_2_adj phq8_severe_2 phq8_severe_2_adj
	using "${tables}/results-mh-mdd-revised-2.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat 2.treat) 
	stats(c_mean c_sd N ctrls pval_ther qval_ther pval_cash qval_cash est_diff pval_diff qval_diff, fmt(3 3 0 0 3 3 3 3) 
	label("Control mean" "Control SD" "Observations" "Controls" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values" "Coeff.: IPT-G+ - IPT-G" "H0: IPT-G=IPT-G+ p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 
* 24m: 
#delimit ;
esttab 	phq8_mdd_3 phq8_mdd_3_adj phq8_severe_3 phq8_severe_3_adj
	using "${tables}/results-mh-mdd-revised-3.tex", replace  booktabs
	label nostar nonote nomtitle nonumber frag nolines nogaps
	b(3) se(3) keep(1.treat 2.treat) 
	stats(c_mean c_sd N ctrls pval_ther qval_ther pval_cash qval_cash est_diff pval_diff qval_diff, fmt(3 3 0 0 3 3 3 3) 
	label("Control mean" "Control SD" "Observations" "Controls" "H0: IPT-G=0 p-values" "\hspace{20pt}FDR adj. q-values" "H0: IPT-G+=0 p-values" "\hspace{20pt}FDR adj. q-values" "Coeff.: IPT-G+ - IPT-G" "H0: IPT-G=IPT-G+ p-values" "\hspace{20pt}FDR adj. q-values")) 
	; 
#delimit cr 

